### Replication files for
### "Place Attachments: Theory and Measurement for Political Science"
### Elizabeth Mitchell Elder and Hans Lueders
### Political Behavior
### July 2025




# Set-up ------------------------------------------------------------------

### packages
library(ggplot2)
library(foreign)
require(lmtest)
library(stargazer)
library(xtable)
library(maptools)
library(gridExtra)
library(ggmap)
library(rgdal)
library(sp)
library(gridExtra)
library(grid)
library(lattice)
library(Matching)
library(ggpubr)
library(geosphere)
library(spatstat)
library(raster)
library(lfe)
library(AER)
library(dplyr)
library(corrplot)
library(ggrepel)
library(rbounds)
library(mgcv)
library(mediation)
library(sensemakr)
library(cjoint)
library(emmeans)
library(cregg)
library(mapproj)
library(readxl)
library(ggpubr)
library(Hmisc)
library(lavaan)
library(estimatr)
library(scales)



### empty environment
rm(list=ls())


### set working directory
path <- "" # define path here
setwd(path)






# Load and prepare datasets -----------------------------------------------

### Germany survey
# load data
load(paste(path, "Data/PlaceAttachments_DEnationwide.Rdata", sep=""))

# create local attachments scale
de$attach_belong <- as.numeric(de$attach_belong)
de$attach_part <- as.numeric(de$attach_part)
de$attach_feelgood <- as.numeric(de$attach_feelgood)
de$attach_memories <- as.numeric(de$attach_memories)
de$attach_stayforever <- as.numeric(de$attach_stayforever)

cfa.mod.de <- cfa('loc.attach =~ attach_belong + attach_part + attach_feelgood + attach_memories + attach_stayforever', data = de, missing="ml")
cfa.de <- summary(cfa.mod.de, fit.measures = TRUE, standardized = TRUE)
de$loc.attach <- rescale(lavPredict(cfa.mod.de, newdata=de)[,1], to=c(0,1))

de$loc.attachDICH <- ifelse(de$loc.attach <= median(de$loc.attach, na.rm=T), 0, 1)

de.conjoint <- merge(de.conjoint, de[,c("ResponseId","loc.attach")], by="ResponseId", all.x=T)

de.conjoint$loc.attachDICH <- ifelse(de.conjoint$loc.attach <= median(de.conjoint$loc.attach, na.rm=T), 0, 1)
de.conjoint$loc.attachDICH <- factor(de.conjoint$loc.attachDICH, levels=0:1,
                              labels = c("weak", "strong"))


### US nationwide survey
# load data
load(paste(path, "Data/PlaceAttachments_USnationwide.Rdata", sep=""))

# create local attachments scale
us$attach_belong <- as.numeric(us$attach_belong)
us$attach_part <- as.numeric(us$attach_part)
us$attach_feelgood <- as.numeric(us$attach_feelgood)
us$attach_memories <- as.numeric(us$attach_memories)
us$attach_stayforever <- as.numeric(us$attach_stayforever)

cfa.mod.us <- cfa('loc.attach =~ attach_belong + attach_part + attach_feelgood + attach_memories + attach_stayforever', data = us, missing="ml")
cfa.us <- summary(cfa.mod.us, fit.measures = TRUE, standardized = TRUE)

us$loc.attach <- rescale(lavPredict(cfa.mod.us, newdata=us)[,1], to=c(0,1))

us$loc.attachDICH <- ifelse(us$loc.attach <= median(us$loc.attach, na.rm=T), 0, 1)

us.conjoint <- merge(us.conjoint, us[,c("ResponseId","loc.attach")], by="ResponseId", all.x=T)

us.conjoint$loc.attachDICH <- ifelse(us.conjoint$loc.attach <= median(us.conjoint$loc.attach, na.rm=T), 0, 1)
us.conjoint$loc.attachDICH <- factor(us.conjoint$loc.attachDICH, levels=0:1,
                              labels = c("weak", "strong"))

### US regional survey
# load data
load(paste(path, "Data/PlaceAttachments_USregional.Rdata", sep=""))

# create local attachments scale
cfa.mod.app <- cfa('loc.attach =~ attach_belong + attach_part + attach_feelgood + attach_memories + attach_stayforever', data = app, missing="ml")
cfa.app <- summary(cfa.mod.app, fit.measures = TRUE, standardized = TRUE)

app$loc.attach <- rescale(lavPredict(cfa.mod.app, newdata=app)[,1], to=c(0,1))

app$loc.attachDICH <- ifelse(app$loc.attach <= median(app$loc.attach, na.rm=T), 0, 1)





# Table 1: Component statements of local attachments scale ----------------

### output container
out <- data.frame(matrix(nrow=5,ncol=4))
colnames(out) <- c("question", "meanDE", "meanUS", "meanApp")
out$question <- factor(1:5, levels = 1:5,
                       labels = c("(1) When I am in my town or city, I feel strongly like I belong there.", 
                                  "(2) Living in my town or city is an important part of who I am.", 
                                  "(3) I feel really at home in my town or city.", 
                                  "(4) If I go through my town or city, I find things that remind me of my past.", 
                                  "(5) I would like to stay in my town or city for the rest of my life."))

### get values
# Germany
out[1,2] <- wtd.mean(de$attach_belong-3, weights = de$weight, na.rm=T)
out[2,2] <- wtd.mean(de$attach_part-3, weights = de$weight, na.rm=T)
out[3,2] <- wtd.mean(de$attach_feelgood-3, weights = de$weight, na.rm=T)
out[4,2] <- wtd.mean(de$attach_memories-3, weights = de$weight, na.rm=T)
out[5,2] <- wtd.mean(de$attach_stayforever-3, weights = de$weight, na.rm=T)

# USA
out[1,3] <- wtd.mean(us$attach_belong-3, weights = us$weight, na.rm=T)
out[2,3] <- wtd.mean(us$attach_part-3, weights = us$weight, na.rm=T)
out[3,3] <- wtd.mean(us$attach_feelgood-3, weights = us$weight, na.rm=T)
out[4,3] <- wtd.mean(us$attach_memories-3, weights = us$weight, na.rm=T)
out[5,3] <- wtd.mean(us$attach_stayforever-3, weights = us$weight, na.rm=T)

# Appalachia
out[1,4] <- wtd.mean(app$attach_belong-3, weights = app$weight, na.rm=T)
out[2,4] <- wtd.mean(app$attach_part-3, weights = app$weight, na.rm=T)
out[3,4] <- wtd.mean(app$attach_feelgood-3, weights = app$weight, na.rm=T)
out[4,4] <- wtd.mean(app$attach_memories-3, weights = app$weight, na.rm=T)
out[5,4] <- wtd.mean(app$attach_stayforever-3, weights = app$weight, na.rm=T)


### Print table
print(xtable(out), include.rownames=FALSE)





# Table 2: Factor loadings and measurement statistics ---------------------

### output container
out <- data.frame(matrix(nrow=9, ncol=4))

colnames(out) <- c("variable", "Germany", "United States", "Appalachia")
out$variable <- c("(1) Belong", "(2) Part of self", "(3) Feel at home", "(4) Memories of past", "(5) Stay forever",
                  "Comparative Fit Index", "Tucker-Lews Index", "Stand. Root Mean Sq. Residual", "Cronbach's Alpha")



### get values
# factor loadings
out[1:5,2] <- inspect(cfa.mod.de,what="std")$lambda
out[1:5,3] <- inspect(cfa.mod.us,what="std")$lambda
out[1:5,4] <- inspect(cfa.mod.app,what="std")$lambda

# measurement statistics -- Germany
out[6,2] <- cfa.de$fit["cfi.robust"]
out[7,2] <- cfa.de$fit["tli.robust"]
out[8,2] <- cfa.de$fit["srmr"]
de2 <- subset(de, !is.na(de$weight))
svy <- survey::svydesign(ids=~1, data=de2, weights=~weight)
out[9,2] <- svycralpha(~attach_belong+attach_part+attach_feelgood+attach_memories+attach_stayforever, svy, na.rm=T)
rm(de2)

# measurement statistics -- USA
out[6,3] <- cfa.us$fit["cfi.robust"]
out[7,3] <- cfa.us$fit["tli.robust"]
out[8,3] <- cfa.us$fit["srmr"]
us2 <- subset(us, !is.na(us$weight))
svy <- survey::svydesign(ids=~1, data=us2, weights=~weight)
out[9,3] <- svycralpha(~attach_belong+attach_part+attach_feelgood+attach_memories+attach_stayforever, svy, na.rm=T)
rm(us2)

# measurement statistics -- Appalachia
out[6,4] <- cfa.app$fit["cfi.robust"]
out[7,4] <- cfa.app$fit["tli.robust"]
out[8,4] <- cfa.app$fit["srmr"]
svy <- survey::svydesign(ids=~1, data=app, weights=~weight)
out[9,4] <- svycralpha(~attach_belong+attach_part+attach_feelgood+attach_memories+attach_stayforever, svy, na.rm=T)


### Table output
print(xtable(out, caption="Factor loadings and measurement statistics", digits=3, align=c("rrccc")), 
      include.rownames=FALSE, caption.placement = "top")





# Figure 1: Distribution of strength of local attachments -----------------

### compute number of statements that a respondent agrees with somewhat or completely
# Germany
de$attach_local_count <- 0
de$attach_local_count[de$attach_belong >= 4 & !is.na(de$attach_belong)] <- de$attach_local_count[de$attach_belong >= 4 & !is.na(de$attach_belong)] + 1
de$attach_local_count[de$attach_feelgood >= 4 & !is.na(de$attach_feelgood)] <- de$attach_local_count[de$attach_feelgood >= 4 & !is.na(de$attach_feelgood)] + 1
de$attach_local_count[de$attach_memories >= 4 & !is.na(de$attach_memories)] <- de$attach_local_count[de$attach_memories >= 4 & !is.na(de$attach_memories)] + 1
de$attach_local_count[de$attach_part >= 4 & !is.na(de$attach_part)] <- de$attach_local_count[de$attach_part >= 4 & !is.na(de$attach_part)] + 1
de$attach_local_count[de$attach_stayforever >= 4 & !is.na(de$attach_stayforever)] <- de$attach_local_count[de$attach_stayforever >= 4 & !is.na(de$attach_stayforever)] + 1

# United States
us$attach_local_count <- 0
us$attach_local_count[us$attach_belong >= 4 & !is.na(us$attach_belong)] <- us$attach_local_count[us$attach_belong >= 4 & !is.na(us$attach_belong)] + 1
us$attach_local_count[us$attach_feelgood >= 4 & !is.na(us$attach_feelgood)] <- us$attach_local_count[us$attach_feelgood >= 4 & !is.na(us$attach_feelgood)] + 1
us$attach_local_count[us$attach_memories >= 4 & !is.na(us$attach_memories)] <- us$attach_local_count[us$attach_memories >= 4 & !is.na(us$attach_memories)] + 1
us$attach_local_count[us$attach_part >= 4 & !is.na(us$attach_part)] <- us$attach_local_count[us$attach_part >= 4 & !is.na(us$attach_part)] + 1
us$attach_local_count[us$attach_stayforever >= 4 & !is.na(us$attach_stayforever)] <- us$attach_local_count[us$attach_stayforever >= 4 & !is.na(us$attach_stayforever)] + 1

# Appalachia
app$attach_local_count <- 0
app$attach_local_count[app$attach_belong >= 4 & !is.na(app$attach_belong)] <- app$attach_local_count[app$attach_belong >= 4 & !is.na(app$attach_belong)] + 1
app$attach_local_count[app$attach_feelgood >= 4 & !is.na(app$attach_feelgood)] <- app$attach_local_count[app$attach_feelgood >= 4 & !is.na(app$attach_feelgood)] + 1
app$attach_local_count[app$attach_memories >= 4 & !is.na(app$attach_memories)] <- app$attach_local_count[app$attach_memories >= 4 & !is.na(app$attach_memories)] + 1
app$attach_local_count[app$attach_part >= 4 & !is.na(app$attach_part)] <- app$attach_local_count[app$attach_part >= 4 & !is.na(app$attach_part)] + 1
app$attach_local_count[app$attach_stayforever >= 4 & !is.na(app$attach_stayforever)] <- app$attach_local_count[app$attach_stayforever >= 4 & !is.na(app$attach_stayforever)] + 1


### summarize
# output container
out <- data.frame(matrix(nrow = 18, ncol=3))
colnames(out) <- c("country", "number", "share")
out$country <- factor(rep(1:3,each=6), levels = 1:3,
                      labels = c("Germany", "United States", "Appalachia"))
out$number <- rep(0:5, 3)

# get values
out[1:6,3] <- prop.table(wtd.table(de$attach_local_count, weights=de$weight)$sum.of.weights)
out[7:12,3] <- prop.table(wtd.table(us$attach_local_count, weights=us$weight)$sum.of.weights)
out[13:18,3] <- prop.table(wtd.table(app$attach_local_count, weights=app$weight)$sum.of.weights)


### plot
p <- ggplot() + 
  geom_bar(data=out, aes(x=number, y=share*100, fill=country),
           stat="identity", position = position_dodge()) +
  annotate("text", x=-0.25, y=3, label = "Germany", size=4.5, angle=90, hjust=0, col="white") +    
  annotate("text", x=0, y=3, label = "United States", size=4.5, angle=90, hjust=0, col="white") +    
  annotate("text", x=0.25, y=3, label = "Appalachia", size=4.5, angle=90, hjust=0, col="white") +      
  xlab("Number of statements") +
  ylab("Share of respondents") +
  scale_x_continuous(breaks=seq(0,5,1)) +
  scale_y_continuous(breaks=seq(0,30,5), limits=c(0,30)) +
  scale_fill_manual(values = c("gray70", "gray40", "gray20"), name="") +  
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title = element_text(size=12, face="bold"),
        legend.position = "none",
        legend.text = element_text(size=12), 
        axis.text = element_text(size=12))     
p
ggsave("Figure1_Attachments_Strength.png",p,width=8,height=5)





# Figure 2: Demographic correlates of local attachments -------------------

### create output containers
out <- data.frame(matrix(nrow=3*44, ncol=5))
colnames(out) <- c("sample", "variable", "level", "b", "se")
out$sample <- factor(rep(1:3, each=44), levels = 1:3,
                     labels = c("Germany", "United States", "Appalachia"))
out$variable <- factor(rep(c(1,1,1,1,1,
                             2,2,
                             3,3,3,3,3,
                             4,4,4,4,
                             5,5,5,5,5,5,
                             6,6,6,6,
                             7,7,7,7,7,
                             8,8,8,8,8,8,
                             9,9,9,9,9,9,9),3), 
                       levels = 1:9,
                       labels = c("Age", "Gender", "Income quintile", 
                                  "Education (DE)", "Education (US)", "Race (US)",
                                  "Time lived in place", "Party preference (DE)", "Party ID (US)"))
out$level <- factor(rep(c(1,2,3,4,5,
                          6,7,
                          8,9,10,11,12,
                          13,14,15,16,
                          17,18,19,20,21,22,
                          23,24,25,26,
                          27,28,29,30,31,
                          32,33,34,35,36,37,
                          38,39,40,41,42,43,44),3),
                    levels = 1:44,
                    labels = c("18-24","25-34", "35-49", "50-64", "65+",
                               "Male", "Female",
                               "1st", "2nd", "3rd", "4th", "5th",
                               "8th grade", "10th grade", "12th/13th grade", "College degree or more",
                               "Less than high school", "High school", "Some college", "2-year college", "4-year college", "Postgraduate degree",
                               "White", "Black", "Latino", "Asian",
                               "0-5 years", "6-10 years", "11-15 years", "16+ years", "Whole life",
                               "Left", "SPD", "Greens", "FDP", "CDU/CSU", "AfD",
                               "Strong Democrat", "Democrat", "Lean Democrat", "Independent", "Lean Republican", "Republican", "Strong Republican"))


### get values
# Germany
out[1:5,4:5] <- summary(lm_robust(loc.attach ~ demo_agegroup - 1, data=de, weights=de$weight))$coef[1:5,1:2]
out[6:7,4:5] <- summary(lm_robust(loc.attach ~ I(demo_female == 1) - 1, data=de, weights=de$weight))$coef[1:2,1:2]
out[8:12,4:5] <- summary(lm_robust(loc.attach ~ demo_income_quintile - 1, data=de, weights=de$weight))$coef[1:5,1:2]
out[13:16,4:5] <- summary(lm_robust(loc.attach ~ demo_educSchoolProf - 1, data=de, weights=de$weight))$coef[1:4,1:2]
out[17:22,4:5] <- NA
out[23:26,4:5] <- NA
out[27:31,4:5] <- summary(lm_robust(loc.attach ~ timelived_place_grouped - 1, data=de, weights=de$weight))$coef[1:5,1:2]
out[32:37,4:5] <- summary(lm_robust(loc.attach ~ party - 1, data=de, weights=de$weight))$coef[1:6,1:2]
out[38:44,4:5] <- NA

# United States
out[45:49,4:5] <- summary(lm_robust(loc.attach ~ demo_agegroup - 1, data=us, weights=us$weight))$coef[1:5,1:2]
out[50:51,4:5] <- summary(lm_robust(loc.attach ~ I(demo_female == 1) - 1, data=us, weights=us$weight))$coef[1:2,1:2]
out[52:56,4:5] <- summary(lm_robust(loc.attach ~ demo_income_quintile - 1, data=us, weights=us$weight))$coef[1:5,1:2]
out[57:60,4:5] <- NA
out[61:66,4:5] <- summary(lm_robust(loc.attach ~ demo_educ - 1, data=us, weights=us$weight))$coef[1:6,1:2]
out[67:70,4:5] <- summary(lm_robust(loc.attach ~ demo_race_recode - 1, data=us, weights=us$weight))$coef[1:4,1:2]
out[71:75,4:5] <- summary(lm_robust(loc.attach ~ timelived_place_grouped - 1, data=us, weights=us$weight))$coef[1:5,1:2]
out[76:81,4:5] <- NA
out[82:88,4:5] <- summary(lm_robust(loc.attach ~ partisanship - 1, data=us, weights=us$weight))$coef[1:7,1:2]

# Appalachia
out[89:93,4:5] <- summary(lm_robust(loc.attach ~ demo_agegroup - 1, data=app, weights=app$weight))$coef[1:5,1:2]
out[94:95,4:5] <- summary(lm_robust(loc.attach ~ I(gender == "F") - 1, data=app, weights=app$weight))$coef[1:2,1:2]
out[96:100,4:5] <- NA
out[101:104,4:5] <- NA
out[105:110,4:5] <- summary(lm_robust(loc.attach ~ as.factor(educ) - 1, data=app, weights=app$weight))$coef[1:6,1:2]
out[111:114,4:5] <- summary(lm_robust(loc.attach ~ demo_race_recode - 1, data=app, weights=app$weight))$coef[1:4,1:2]
out[115:119,4:5] <- NA
out[120:125,4:5] <- NA
out[126:132,4:5] <- summary(lm_robust(loc.attach ~ partisanship - 1, data=app, weights=app$weight))$coef[1:7,1:2]


### Plot
p <- ggplot() + 
  geom_point(data=out,aes(x=level,y=b, col=sample, shape=sample),
             size=2.5, position=position_dodge(width=0.6)) +
  geom_errorbar(data=out, aes(x=level, ymin=b-1.96*se, ymax=b+1.96*se, col=sample), width=0.25,
                linewidth=1, position=position_dodge(width=0.6)) +
  scale_shape_manual(values = c(15,17,19), name="") +
  scale_color_manual(values = c("gray40","burlywood3", "burlywood4"), name="") +
  ylab("Strength of Local Attachment") +
  coord_flip() +
  facet_wrap(~variable, scales="free_y") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size=12), 
        axis.text = element_text(size=12))   
p

ggsave("Figure2_Attachments_Demographics.png",p,width=12,height=7)






# Figure 3: Place-based correlates of local attachments -------------------

### create output container
out <- data.frame(matrix(nrow=6*3*3, ncol=5))
colnames(out) <- c("country", "outcome", "level", "b", "se")
out$country <- factor(rep(1:3,each=18), levels = 1:3,
                      labels = c("Germany", "United States", "Appalachia"))
out$outcome <- factor(rep(rep(1:6,each=3),3), levels = 1:6,
                      labels = c("Rural/urban", "Income","Unemployment",
                                 "Life expectancy","Homicide rate",
                                 "Deaths of despair"))
out$level <- factor(rep(c(1:3,rep(4:6,5)),3), levels = 1:6,
                    labels = c("Rural/non-metro", "Urban/metro", "NA", "Low", "Medium", "High"))


### get values
# Germany
out[1:2,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_urban) - 1, data=de, weights=de$weight))$coef[1:2,1:2]
out[4:6,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_hhincome_bin) - 1, data=de, weights=de$weight))$coef[1:3,1:2]
out[7:9,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_uer_bin) - 1, data=de, weights=de$weight))$coef[1:3,1:2]

# United States
out[19:20,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_metro) - 1, data=us, weights=us$weight))$coef[1:2,1:2]
out[22:24,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_hhincome_bin) - 1, data=us, weights=us$weight))$coef[1:3,1:2]
out[25:27,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_uer_bin) - 1, data=us, weights=us$weight))$coef[1:3,1:2]
out[28:30,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_le14_bin) - 1, data=us, weights=us$weight))$coef[1:3,1:2]
out[31:33,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_homicides_bin) - 1, data=us, weights=us$weight))$coef[1:3,1:2]

# Appalachia
out[37:38,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_metro) - 1, data=app, weights=app$weight))$coef[1:2,1:2]
out[40:42,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_hhincome_bin) - 1, data=app, weights=app$weight))$coef[1:3,1:2]
out[43:45,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_uer_bin) - 1, data=app, weights=app$weight))$coef[1:3,1:2]
out[46:48,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_le14_bin) - 1, data=app, weights=app$weight))$coef[1:3,1:2]
out[49:51,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_homicides_bin) - 1, data=app, weights=app$weight))$coef[1:3,1:2]
out[52:54,4:5] <- summary(lm_robust(loc.attach ~ as.factor(county_deaths_all_bin) - 1, data=app, weights=app$weight))$coef[1:3,1:2]


### Plot
p <- ggplot() + 
  geom_point(data=out[out$level != "NA",], aes(x=level, y=b, shape=country, col=country),
             position=position_dodge(width=0.5), size=2.5) +
  geom_errorbar(data=out[out$level != "NA",], aes(x=level, ymin=b-1.96*se, ymax=b+1.96*se, col=country),
                position=position_dodge(width=0.5), linewidth=1, width=0.25) +
  facet_wrap(~outcome, scales="free_y") +
  scale_shape_manual(values = c(15,17,19), name="") +
  scale_color_manual(values = c("gray40","burlywood3", "burlywood4"), name="") +
  ylab("Strength of Local Attachment") +
  coord_flip() +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size=12, face="bold"),
        legend.position = "bottom",
        legend.text = element_text(size=12), 
        axis.text = element_text(size=12))   
p
ggsave("Figure3_Attachments_Econsit.png", p, width=10, height=5)





# Figure 4: Local attachments are distinct from rural and urban id --------


### Correlation with rural/urban identity
# create output container
means <- data.frame(matrix(nrow=10,ncol=5))
colnames(means) <- c("urban","identity_rururb", "mean", "sd", "n")
means$urban <- rep(0:1, each=5)
means$urban <- factor(means$urban, levels = 0:1,
                      labels = c("Rural", "Urban"))
means$identity_rururb <- factor(1:5, levels = 1:5,
                                labels = c("none at all", "a little", "a moderate amount",
                                           "a lot", "a great deal"))

# get means
means$mean[1] <- wtd.mean(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "none at all"], 
                          weights = us$weight[us$demo_urban == "Rural" & us$identity_rururb == "none at all"],
                          na.rm=T)
means$mean[2] <- wtd.mean(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a little"], 
                          weights = us$weight[us$demo_urban == "Rural" & us$identity_rururb == "a little"],
                          na.rm=T)
means$mean[3] <- wtd.mean(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a moderate amount"], 
                          weights = us$weight[us$demo_urban == "Rural" & us$identity_rururb == "a moderate amount"],
                          na.rm=T)
means$mean[4] <- wtd.mean(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a lot"], 
                          weights = us$weight[us$demo_urban == "Rural" & us$identity_rururb == "a lot"],
                          na.rm=T)
means$mean[5] <- wtd.mean(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a great deal"], 
                          weights = us$weight[us$demo_urban == "Rural" & us$identity_rururb == "a great deal"],
                          na.rm=T)

means$mean[6] <- wtd.mean(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "none at all"], 
                          weights = us$weight[us$demo_urban == "Urban" & us$identity_rururb == "none at all"],
                          na.rm=T)
means$mean[7] <- wtd.mean(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a little"], 
                          weights = us$weight[us$demo_urban == "Urban" & us$identity_rururb == "a little"],
                          na.rm=T)
means$mean[8] <- wtd.mean(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a moderate amount"], 
                          weights = us$weight[us$demo_urban == "Urban" & us$identity_rururb == "a moderate amount"],
                          na.rm=T)
means$mean[9] <- wtd.mean(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a lot"], 
                          weights = us$weight[us$demo_urban == "Urban" & us$identity_rururb == "a lot"],
                          na.rm=T)
means$mean[10] <- wtd.mean(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a great deal"], 
                           weights = us$weight[us$demo_urban == "Urban" & us$identity_rururb == "a great deal"],
                           na.rm=T)

# get SD
means$sd[1] <- wtd.var(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "none at all"], 
                       weights = us$weight[us$demo_urban == "Rural" & us$identity_rururb == "none at all"],
                       na.rm=T)
means$sd[2] <- wtd.var(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a little"], 
                       weights = us$weight[us$demo_urban == "Rural" & us$identity_rururb == "a little"],
                       na.rm=T)
means$sd[3] <- wtd.var(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a moderate amount"], 
                       weights = us$weight[us$demo_urban == "Rural" & us$identity_rururb == "a moderate amount"],
                       na.rm=T)
means$sd[4] <- wtd.var(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a lot"], 
                       weights = us$weight[us$demo_urban == "Rural" & us$identity_rururb == "a lot"],
                       na.rm=T)
means$sd[5] <- wtd.var(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a great deal"], 
                       weights = us$weight[us$demo_urban == "Rural" & us$identity_rururb == "a great deal"],
                       na.rm=T)

means$sd[6] <- wtd.var(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "none at all"], 
                       weights = us$weight[us$demo_urban == "Urban" & us$identity_rururb == "none at all"],
                       na.rm=T)
means$sd[7] <- wtd.var(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a little"], 
                       weights = us$weight[us$demo_urban == "Urban" & us$identity_rururb == "a little"],
                       na.rm=T)
means$sd[8] <- wtd.var(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a moderate amount"], 
                       weights = us$weight[us$demo_urban == "Urban" & us$identity_rururb == "a moderate amount"],
                       na.rm=T)
means$sd[9] <- wtd.var(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a lot"], 
                       weights = us$weight[us$demo_urban == "Urban" & us$identity_rururb == "a lot"],
                       na.rm=T)
means$sd[10] <- wtd.var(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a great deal"], 
                        weights = us$weight[us$demo_urban == "Urban" & us$identity_rururb == "a great deal"],
                        na.rm=T)
means$sd <- sqrt(means$sd)

# number of observations
means$n[1] <- length(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "none at all"])
means$n[2] <- length(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a little"])
means$n[3] <- length(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a moderate amount"])
means$n[4] <- length(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a lot"])
means$n[5] <- length(us$loc.attach[us$demo_urban == "Rural" & us$identity_rururb == "a great deal"])

means$n[6] <- length(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "none at all"])
means$n[7] <- length(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a little"])
means$n[8] <- length(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a moderate amount"])
means$n[9] <- length(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a lot"])
means$n[10] <- length(us$loc.attach[us$demo_urban == "Urban" & us$identity_rururb == "a great deal"])

# compute standard error
means$se <- means$sd / sqrt(means$n)


## plot
p <- ggplot() + 
  geom_jitter(data=us[!is.na(us$loc.attach) & !is.na(us$identity_rururb),], 
              aes(x=as.numeric(identity_rururb), y=loc.attach),
              col="gray20", fill="gray20", alpha=.65) +
  geom_line(data=means, aes(x=as.numeric(identity_rururb), y=mean), linewidth=1.5, col="chocolate4") +
  geom_ribbon(data=means, aes(x=as.numeric(identity_rururb), ymin=mean-1.96*se, ymax=mean+1.96*se),
              alpha=.7, fill="burlywood2") +
  scale_x_continuous(breaks=1:5,
                     labels = c("none\nat all", "a little", "a moderate\namount", "a lot", "a great\ndeal")) +
  xlab("Strength of Rural/Urban Identity") +
  ylab("Strength of Local Attachment") +
  facet_wrap(~urban) +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title = element_text(size=12, face="bold"),
        legend.position = "bottom",
        legend.text = element_text(size=12), 
        axis.text = element_text(size=12))     
p
ggsave("Figure4_Attachments_RurUrbIdentityUS.png",p,width=9,height=5)





# Figure 5: Correlations between local attachments and engagement  --------

### Write a function
regfun <- function(dv){
  # define outcome
  de$dv <- as.numeric(de[,dv])
  us$dv <- as.numeric(us[,dv])
  
  # regressions
  outDE <- summary(lm_robust(dv ~ loc.attach, data=de, weights=de$weight))$coef[2,1:2]
  outUS <- summary(lm_robust(dv ~ loc.attach, data=us, weights=us$weight))$coef[2,1:2]
  
  # convert to SD changes
  outDE.SD <- sqrt(weighted.var(de$loc.attach, na.rm=T, w=de$weight))*outDE / sqrt(weighted.var(de$dv, na.rm=T, w=de$weight))
  outUS.SD <- sqrt(weighted.var(us$loc.attach, na.rm=T, w=us$weight))*outUS / sqrt(weighted.var(us$dv, na.rm=T, w=us$weight))
  
  # return
  return(rbind(outDE.SD, outUS.SD))
}


### output container
out <- data.frame(matrix(nrow=42,ncol=6))
colnames(out) <- c("analysis","outcome", "level", "country", "b", "se")
out$analysis <- factor(rep(c(1,1,1,1,2,2), 7), levels = 1:2,
                       labels = c("level", "Difference federal vs. local"))
out$outcome <- factor(rep(7:1, each=6), levels = 1:7,
                      labels = c("Identification",
                                 "Importance of turnout",
                                 "News consumption",
                                 "Discuss with family",
                                 "Gov't has impact on respondent",
                                 "Respondent can influence gov't",
                                 "Trust in gov't"))
out$level <- factor(rep(rep(1:3,each=2), 7), levels = 1:3,
                    labels = c("Federal", "Local", "Difference"))
out$country <- factor(rep(1:2), levels = 1:2,
                      labels = c("Germany", "United States"))


### Get values
out[1:2,5:6] <- regfun(dv = "govlev_trust_federal")
out[3:4,5:6] <- regfun(dv = "govlev_trust_local")
out[5:6,5:6] <- regfun(dv = "govlev_trustDiff_local")

out[7:8,5:6] <- regfun(dv = "efficacy_federal")
out[9:10,5:6] <- regfun(dv = "efficacy_local")
out[11:12,5:6] <- regfun(dv = "efficacyDiff_local")

out[13:14,5:6] <- regfun(dv = "govlev_effect_federal")
out[15:16,5:6] <- regfun(dv = "govlev_effect_local")
out[17:18,5:6] <- regfun(dv = "govlev_effectDiff_local")

out[19:20,5:6] <- regfun(dv = "discuss_federal")
out[21:22,5:6] <- regfun(dv = "discuss_local")
out[23:24,5:6] <- regfun(dv = "discuss_difference")

out[25:26,5:6] <- regfun(dv = "news_federal")
out[27:28,5:6] <- regfun(dv = "news_local")
out[29:30,5:6] <- regfun(dv = "news_difference")

out[31:32,5:6] <- NA
out[33:34,5:6] <- NA
out[35:36,5:6] <- regfun(dv = "turnout_level")

out[37:38,5:6] <- NA
out[39:40,5:6] <- NA
out[41:42,5:6] <- regfun(dv = "identity_country_local")



### Plot
p1 <- ggplot() +
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out[out$analysis=="level",], aes(x=outcome, y=b, col=level, shape=level),
             position = position_dodge(width=0.6), size=2.5) +
  geom_errorbar(data=out[out$analysis=="level",], aes(x=outcome, ymin=b-1.96*se, ymax=b+1.96*se, col=level),
                position = position_dodge(width=0.6), linewidth=1, width=0.25) +
  scale_y_continuous(limits = c(-0.1,0.4)) +  
  facet_wrap(~country) +
  scale_color_manual(values=c("burlywood3", "burlywood4"), name="") +
  scale_shape_manual(values=c(19,17), name="") +
  ylab("Estimated relationship with local attachments\n(standard deviations)") +
  coord_flip() +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size=12), 
        axis.text = element_text(size=12))  

p2 <- ggplot() +
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out[out$analysis=="Difference federal vs. local",], aes(x=outcome, y=b, col=country, shape=country),
             position = position_dodge(width=0.6), size=2.5) +
  geom_errorbar(data=out[out$analysis=="Difference federal vs. local",], aes(x=outcome, ymin=b-1.96*se, ymax=b+1.96*se, col=country),
                position = position_dodge(width=0.6), linewidth=1, width=0.25) +
  scale_y_continuous(limits = c(-0.3,0.06)) +
  facet_wrap(~analysis) +
  scale_color_manual(values=c("chocolate2", "chocolate4"), name="") +
  scale_shape_manual(values=c(19,17), name="") +
  ylab("Estimated relationship with local attachments\n(standard deviations)") +
  coord_flip() +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size=12), 
        axis.text.x = element_text(size=12),
        axis.text.y = element_blank())  

p <- ggarrange(p1,p2,widths=c(0.72,0.28))
p

ggsave("Figure5_Attachments_Differences.png",p,width=10,height=6)






# Figure 6: Local considerations have a stronger effect among vote --------

### prepare output container
out <- data.frame(matrix(nrow=88, ncol=7))
colnames(out) <- c("country", "sample", "dimension", "attribute", "b", "lb", "ub")
out$country <- factor(rep(1:4,each=22), levels = 1:4,
                      labels = c("Germany", "United States", "Difference in Marginal Means", "Difference in AMCEs"))
out$sample <- factor(c(rep(rep(1:2,each=11),2), rep(3:4, each=11), rep(3:4, each=11)), levels = 1:4,
                     labels = c("Strong", "Weak", "Germany", "United States"))
out$dimension <- factor(rep(c(1,1,1,2,2,2,2,3,3,3,3),8), levels = 1:3,
                        labels = c("Position on federalism", "Experience in local politics", "Time lived in place"))
out$attribute <- factor(rep(11:1,8), levels = 1:11,
                        labels = c("20+ years   ", "10-15 years   ", "1-5 years   ", "Time lived in place",
                                   "9-15 years   ","1-5  years   ", "None   ", "Experience in local politics", 
                                   "Rather local   ","Rather federal   ", "Position on federalism"))


### Marginal means
# Germany
reg <- cj(data = de.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + C_party + 
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "mm",
          by=~loc.attachDICH)

out[2:3,5:7] <- reg[43:44,c(6,10,11)]
out[5:7,5:7] <- reg[40:42,c(6,10,11)]
out[9:11,5:7] <- reg[37:39,c(6,10,11)]

out[13:14,5:7] <- reg[21:22,c(6,10,11)]
out[16:18,5:7] <- reg[18:20,c(6,10,11)]
out[20:22,5:7] <- reg[15:17,c(6,10,11)]


# United States
reg <- cj(data = us.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + 
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "mm",
          by=~loc.attachDICH)

out[24:25,5:7] <- reg[37:38,c(6,10,11)]
out[27:29,5:7] <- reg[34:36,c(6,10,11)]
out[31:33,5:7] <- reg[31:33,c(6,10,11)]

out[35:36,5:7] <- reg[18:19,c(6,10,11)]
out[38:40,5:7] <- reg[15:17,c(6,10,11)]
out[42:44,5:7] <- reg[12:14,c(6,10,11)]



### Difference in marginal means
# Germany
reg <- cj(data = de.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + C_party +
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "mm_difference",
          by=~loc.attachDICH)

out[46:47,5:7] <- reg[21:22,c(6,10,11)]
out[49:51,5:7] <- reg[18:20,c(6,10,11)]
out[53:55,5:7] <- reg[15:17,c(6,10,11)]

# United States
reg <- cj(data = us.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + 
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "mm_difference",
          by=~loc.attachDICH)

out[57:58,5:7] <- reg[18:19,c(6,10,11)]
out[60:62,5:7] <- reg[15:17,c(6,10,11)]
out[64:66,5:7] <- reg[12:14,c(6,10,11)]




### Difference in AMCEs
# Germany
reg <- cj(data = de.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + C_party +
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "amce_difference",
          by=~loc.attachDICH)

out[68,5] <- 0
out[69,5:7] <- reg[4,c(7,11,12)]
out[71,5] <- 0
out[72:73,5:7] <- reg[6:7,c(7,11,12)]
out[75,5] <- 0
out[76:77,5:7] <- reg[14:15,c(7,11,12)]


# United States
reg <- cj(data = us.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + 
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "amce_difference",
          by=~loc.attachDICH)

out[79,5] <- 0
out[80,5:7] <- reg[4,c(7,11,12)]
out[82,5] <- 0
out[83:84,5:7] <- reg[6:7,c(7,11,12)]
out[86,5] <- 0
out[87:88,5:7] <- reg[12:13,c(7,11,12)]



### create new variables for faceting
out$facet_country <- NA
out$facet_country[out$country == "Germany" | out$sample == "Germany"] <- 1
out$facet_country[out$country == "United States" | out$sample == "United States"] <- 2
out$facet_country <- factor(out$facet_country, levels = 1:2,
                            labels = c("Germany", "United States"))

out$facet_analysis <- NA
out$facet_analysis[out$country %in% c("Germany", "United States")] <- 1
out$facet_analysis[out$country == "Difference in Marginal Means"] <- 2
out$facet_analysis[out$country == "Difference in AMCEs"] <- 3
out$facet_analysis <- factor(out$facet_analysis, levels = 1:3,
                             labels = c("Marginal Means", "Difference in\nMarginal Means", "Difference in \nAMCEs"))

out$hline <- ifelse(out$facet_analysis == "Marginal Means", 0.5, 0)

out$sample2 <- as.numeric(out$sample)
out$sample2[out$sample2 == 4] <- 3
out$sample2 <- factor(out$sample2, levels = 1:3,
                      labels = c("Strong", "Weak", "Strong - Weak"))


### Plot
p <- ggplot() + 
  geom_hline(data=out, aes(yintercept = hline), lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, shape=sample2, col=sample2), 
             size=3, position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, col=sample2), width=0,
                position=position_dodge(width=0.8), linewidth=1) + 
  scale_color_manual(values = c("burlywood4","burlywood3", "gray40"), name="") +
  scale_shape_manual(values = c(15,17,19), name="") +
  facet_grid(facet_country ~ facet_analysis, scales="free_x") +
  coord_flip() + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=13,
                                   face = c(rep("plain",3),"bold",
                                            rep("plain",3), "bold",
                                            rep("plain",2), "bold")),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))
p
ggsave("Figure6_Attachments_Conjoint.png",p,width=11,height=6)







# Figure A1: Screeplot ----------------------------------------------------

### load psych package
library(psych)

### output container
out <- data.frame(matrix(nrow=15, ncol=3))
colnames(out) <- c("country", "index", "eigenvalue")
out$country <- factor(rep(1:3,each=5), levels=1:3,
                      labels = c("Germany", "United States", "Appalachia"))
out$index <- rep(1:5,3)

### get eigenvalues
out[1:5,3] <- (fa(app[,c("attach_belong","attach_part","attach_feelgood","attach_memories","attach_stayforever")], nfactors = 5, 
                  rotate = "none", fm = "ml"))$e.values

out[6:10,3] <- (fa(us[,c("attach_belong","attach_part","attach_feelgood","attach_memories","attach_stayforever")], nfactors = 5, 
                   rotate = "none", fm = "ml"))$e.values

out[11:15,3] <- (fa(de[,c("attach_belong","attach_part","attach_feelgood","attach_memories","attach_stayforever")], nfactors = 5, 
                    rotate = "none", fm = "ml"))$e.values


### Plot
p <- ggplot() + 
  geom_hline(yintercept = 1, lty=2) +
  geom_point(data=out, aes(x=index, y=eigenvalue, col=country), size=2.5) +
  geom_line(data=out, aes(x=index, y=eigenvalue, col=country), linewidth=1) +
  scale_color_manual(values = c("gray20", "gray50", "gray75"), name="") +
  scale_y_continuous(breaks = seq(0,3.5,0.5)) +
  facet_wrap(~country) +
  xlab("Factor index") +
  ylab("Eigenvalue") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title = element_text(size=12, face="bold"),
        legend.position = "none",
        legend.text = element_text(size=12), 
        axis.text = element_text(size=12))   
p
ggsave("FigureA1_Screeplot.png", p, width=8, height=3)





# Tables A1-A3: Correlation matrices --------------------------------------

### Germany
out <- data.frame(round(cor(de[,c("attach_belong","attach_part","attach_feelgood","attach_memories","attach_stayforever")], use="complete.obs"),3))
colnames(out) <- c("Belong", "Part of self", "Feel at home", "Memories of past", "Stay forever")
rownames(out) <- c("Belong", "Part of self", "Feel at home", "Memories of past", "Stay forever")
print(xtable(out, caption="Correlation matrix: Germany survey", digits=3, align=c("rccccc")), caption.placement = "top")

### United States
out <- data.frame(round(cor(us[,c("attach_belong","attach_part","attach_feelgood","attach_memories","attach_stayforever")], use="complete.obs"),3))
colnames(out) <- c("Belong", "Part of self", "Feel at home", "Memories of past", "Stay forever")
rownames(out) <- c("Belong", "Part of self", "Feel at home", "Memories of past", "Stay forever")
print(xtable(out, caption="Correlation matrix: United States survey", digits=3, align=c("rccccc")), caption.placement = "top")

### Appalachia
out <- data.frame(round(cor(app[,c("attach_belong","attach_part","attach_feelgood","attach_memories","attach_stayforever")], use="complete.obs"),3))
colnames(out) <- c("Belong", "Part of self", "Feel at home", "Memories of past", "Stay forever")
rownames(out) <- c("Belong", "Part of self", "Feel at home", "Memories of past", "Stay forever")
print(xtable(out, caption="Correlation matrix: Appalachia survey", digits=3, align=c("rccccc")), caption.placement = "top")




# Figure A2: Political correlates of local attachments in Appalach --------

### Write a function
regfun <- function(dv){
  # define outcome
  app$dv <- as.numeric(app[,dv])
  
  # regressions
  out <- summary(lm_robust(dv ~ loc.attach, data=app, weights=app$weight))$coef[2,1:2]
  
  # convert to SD changes
  out.SD <- sqrt(weighted.var(app$loc.attach, na.rm=T, w=app$weight))*out / sqrt(weighted.var(app$dv, na.rm=T, w=app$weight))
  
  # return
  return(out.SD)
}


### output container
out <- data.frame(matrix(nrow=28,ncol=4))
colnames(out) <- c("class", "outcome", "b", "se")
out$class <- factor(c(rep(1,14), rep(2,6), rep(3,3), rep(4,5)),
                    levels = 1:4,
                    labels = c("Trust", "Engagement", "Efficacy", "Populism and Democracy"))
out$outcome <- factor(1:28, levels = 1:28,
                      labels = c("Trust in Federal Gov't",
                                 "Trust in State Gov't",
                                 "Trust in Local Gov't",
                                 "Trust Gov't to be: Honest",
                                 "Trust Gov't to be: Fair",
                                 "Trust Gov't to be: Competent",
                                 "Trust in Bureaucracy",
                                 "Trust in Charity",
                                 "Trust in Church",
                                 "Trust in Courts",
                                 "Trust in Media",
                                 "Trust in Military",
                                 "Trust in Police",
                                 "Trust in Unions",
                                 "Interest in politics",
                                 "Worked with community",
                                 "Voted",
                                 "Volunteered",
                                 "Discuss politics",
                                 "Contacted official",
                                 "Efficacy: Powerful",
                                 "Efficacy: Helpless",
                                 "Efficacy: Collective",
                                 "Populism: Popular Will",
                                 "Populism: Anti-elitism",
                                 "Manicheism",
                                 "Conspiracism",
                                 "Antidemocratic attitudes"))

### Get values
out[1,3:4] <- regfun(dv = "trust_anes_fed")
out[2,3:4] <- regfun(dv = "trust_anes_state")
out[3,3:4] <- regfun(dv = "trust_anes_local")
out[4,3:4] <- regfun(dv = "trust_factor_corr")
out[5,3:4] <- regfun(dv = "trust_factor_fair")
out[6,3:4] <- regfun(dv = "trust_factor_comp")
out[7,3:4] <- regfun(dv = "trust_orgs_bureau")
out[8,3:4] <- regfun(dv = "trust_orgs_charity")
out[9,3:4] <- regfun(dv = "trust_orgs_church")
out[10,3:4] <- regfun(dv = "trust_orgs_courts")
out[11,3:4] <- regfun(dv = "trust_orgs_media")
out[12,3:4] <- regfun(dv = "trust_orgs_armed_forces")
out[13,3:4] <- regfun(dv = "trust_orgs_police")
out[14,3:4] <- regfun(dv = "trust_orgs_unions")
out[15,3:4] <- regfun(dv = "part_interest")
out[16,3:4] <- regfun(dv = "part_worktog")
out[17,3:4] <- regfun(dv = "part_vote")
out[18,3:4] <- regfun(dv = "part_vol")
out[19,3:4] <- regfun(dv = "part_talk")
out[20,3:4] <- regfun(dv = "part_contact")
out[21,3:4] <- regfun(dv = "powerful")
out[22,3:4] <- regfun(dv = "helpless")
out[23,3:4] <- regfun(dv = "collective")
out[24,3:4] <- regfun(dv = "popwill")
out[25,3:4] <- regfun(dv = "antielite")
out[26,3:4] <- regfun(dv = "manich1")
out[27,3:4] <- regfun(dv = "conspiracy1")
out[28,3:4] <- regfun(dv = "antidem")



### Plot
p <- ggplot() +
  annotate("rect", ymin=-Inf, ymax=Inf,xmin=14.5,xmax=28.5, fill="burlywood2", alpha=.5) +
  annotate("rect", ymin=-Inf, ymax=Inf,xmin=8.5,xmax=14.5, fill="burlywood3", alpha=.5) +
  annotate("rect", ymin=-Inf, ymax=Inf,xmin=5.5,xmax=8.5, fill="burlywood4", alpha=.5) +
  annotate("rect", ymin=-Inf, ymax=Inf,xmin=0.5,xmax=5.5, fill="gray40", alpha=.5) +
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out, aes(x=outcome, y=b),
             position = position_dodge(width=0.6), size=2.5) +
  geom_errorbar(data=out, aes(x=outcome, ymin=b-1.96*se, ymax=b+1.96*se),
                position = position_dodge(width=0.6), linewidth=1, width=0.25) +
  scale_x_discrete(limits=rev) +
  scale_color_manual(values=c("gray30", "gray70"), name="") +
  scale_shape_manual(values=c(19,17), name="") +
  ylab("Estimated relationship with local attachments\n(standard deviations)") +
  coord_flip() +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size=12), 
        axis.text = element_text(size=12))
p

ggsave("FigureA2_Attachments_Appalachia.png", p, width=8,height=8)







# Figure A3: Correlations with basic controls -----------------------------

### Write a function
regfun <- function(dv){
  # define outcome
  de$dv <- as.numeric(de[,dv])
  us$dv <- as.numeric(us[,dv])
  
  # regressions
  outDE <- summary(lm_robust(dv ~ loc.attach + demo_agegroup + demo_female + demo_income_quintile + demo_educSchoolProf + timelived_place_grouped, 
                             data=de, weights=de$weight))$coef[2,1:2]
  outUS <- summary(lm_robust(dv ~ loc.attach + demo_agegroup + demo_female + demo_income_quintile + demo_educ + demo_race_recode + timelived_place_grouped, 
                             data=us, weights=us$weight))$coef[2,1:2]
  
  # convert to SD changes
  outDE.SD <- sqrt(weighted.var(de$loc.attach, na.rm=T, w=de$weight))*outDE / sqrt(weighted.var(de$dv, na.rm=T, w=de$weight))
  outUS.SD <- sqrt(weighted.var(us$loc.attach, na.rm=T, w=us$weight))*outUS / sqrt(weighted.var(us$dv, na.rm=T, w=us$weight))
  
  # return
  return(rbind(outDE.SD, outUS.SD))
}


### output container
out <- data.frame(matrix(nrow=42,ncol=6))
colnames(out) <- c("analysis","outcome", "level", "country", "b", "se")
out$analysis <- factor(rep(c(1,1,1,1,2,2), 7), levels = 1:2,
                       labels = c("level", "Difference federal vs. local"))
out$outcome <- factor(rep(7:1, each=6), levels = 1:7,
                      labels = c("Identification",
                                 "Importance of turnout",
                                 "News consumption",
                                 "Discuss with family",
                                 "Gov't has impact on respondent",
                                 "Respondent can influence gov't",
                                 "Trust in gov't"))
out$level <- factor(rep(rep(1:3,each=2), 7), levels = 1:3,
                    labels = c("Federal", "Local", "Difference"))
out$country <- factor(rep(1:2), levels = 1:2,
                      labels = c("Germany", "United States"))


### Get values
out[1:2,5:6] <- regfun(dv = "govlev_trust_federal")
out[3:4,5:6] <- regfun(dv = "govlev_trust_local")
out[5:6,5:6] <- regfun(dv = "govlev_trustDiff_local")

out[7:8,5:6] <- regfun(dv = "efficacy_federal")
out[9:10,5:6] <- regfun(dv = "efficacy_local")
out[11:12,5:6] <- regfun(dv = "efficacyDiff_local")

out[13:14,5:6] <- regfun(dv = "govlev_effect_federal")
out[15:16,5:6] <- regfun(dv = "govlev_effect_local")
out[17:18,5:6] <- regfun(dv = "govlev_effectDiff_local")

out[19:20,5:6] <- regfun(dv = "discuss_federal")
out[21:22,5:6] <- regfun(dv = "discuss_local")
out[23:24,5:6] <- regfun(dv = "discuss_difference")

out[25:26,5:6] <- regfun(dv = "news_federal")
out[27:28,5:6] <- regfun(dv = "news_local")
out[29:30,5:6] <- regfun(dv = "news_difference")

out[31:32,5:6] <- NA
out[33:34,5:6] <- NA
out[35:36,5:6] <- regfun(dv = "turnout_level")

out[37:38,5:6] <- NA
out[39:40,5:6] <- NA
out[41:42,5:6] <- regfun(dv = "identity_country_local")



### Plot
p1 <- ggplot() +
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out[out$analysis=="level",], aes(x=outcome, y=b, col=level, shape=level),
             position = position_dodge(width=0.6), size=2.5) +
  geom_errorbar(data=out[out$analysis=="level",], aes(x=outcome, ymin=b-1.96*se, ymax=b+1.96*se, col=level),
                position = position_dodge(width=0.6), linewidth=1, width=0.25) +
  scale_y_continuous(limits = c(-0.15,0.4), breaks = seq(-0.1,0.4,0.1)) +  
  facet_wrap(~country) +
  scale_color_manual(values=c("burlywood3", "burlywood4"), name="") +
  scale_shape_manual(values=c(19,17), name="") +
  ylab("Estimated relationship with local attachments\n(standard deviations)") +
  coord_flip() +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        #axis.title.x = element_text(size=14, face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size=12), 
        axis.text = element_text(size=12))  

p2 <- ggplot() +
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out[out$analysis=="Difference federal vs. local",], aes(x=outcome, y=b, col=country, shape=country),
             position = position_dodge(width=0.6), size=2.5) +
  geom_errorbar(data=out[out$analysis=="Difference federal vs. local",], 
                aes(x=outcome, ymin=b-1.96*se, ymax=b+1.96*se, col=country),
                position = position_dodge(width=0.6), linewidth=1, width=0.25) +
  scale_y_continuous(limits = c(-0.3,0.06)) +
  facet_wrap(~analysis) +
  scale_color_manual(values=c("chocolate2", "chocolate4"), name="") +
  scale_shape_manual(values=c(19,17), name="") +
  ylab("Estimated relationship with local attachments\n(standard deviations)") +
  coord_flip() +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        #axis.title.x = element_text(size=14, face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size=12), 
        axis.text.x = element_text(size=12),
        axis.text.y = element_blank())  

p <- ggarrange(p1,p2,widths=c(0.72,0.28))
p

ggsave("FigureA3_Attachments_Differences_democontrols.png",p,width=10,height=6)








# Figure A4: Conjoint full results ----------------------------------------

### full conjoint results
# regressions
regDE <- cj(data = de.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + C_party +
              C_timelived + C_localexp + C_federalism, id = ~ResponseId,
            estimate = "mm",
            by=~loc.attachDICH)
regUS <- cj(data = us.conjoint, formula = C_choice ~ C_age + C_gender + C_occup +
              C_timelived + C_localexp + C_federalism, id = ~ResponseId,
            estimate = "mm",
            by=~loc.attachDICH)

# prepare output container
out <- data.frame(matrix(nrow=136, ncol=7))
colnames(out) <- c("country", "sample", "dimension", "attribute", "b", "lb", "ub")
out$country <- factor(c(rep(1,72), rep(2,64)), levels = 1:2,
                      labels = c("Germany", "United States"))
out$sample <- factor(c(rep(1:2,each=36), rep(1:2,each=32)), levels = 1:2,
                     labels = c("Weak", "Strong"))
out$dimension <- factor(c(rep(c(1,1,1,1,1,
                                2,2,2,2,2,2,
                                3,3,3,3,3,3,3,
                                4,4,4,4,
                                5,5,5,5,5,5,
                                6,6,6,
                                7,7,7,7,7),2),
                          rep(c(1,1,1,1,1,
                                2,2,2,2,2,2,
                                3,3,3,3,3,3,3,
                                5,5,5,5,5,5,
                                6,6,6,
                                7,7,7,7,7),2)), 
                        levels = 1:7,
                        labels = c("Position on federalism", "Experience in local politics", "Time lived in place",
                                   "Party affiliation", "Occupation", "Gender", "Age"))
out$attribute <- c(rep(36:1,2),
                   rep(c(36:19,14:1),2))
out$attribute <- factor(out$attribute, levels = 1:36,
                        labels = c("59 years   ", "49 years   ", "39 years   ", "29 years   ", "Age",
                                   "Female   ", "Male   ", "Gender",
                                   "Attorney   ", "Police officer   ", "Small business owner   ", "Teacher   ", "Civil servant   ", "Occupation",
                                   "Independent local party   ", "SPD   ", "CDU/CSU   ", "Party affilication",
                                   "Whole life   ","20+ years   ", "15 years   ", "10 years   ", "5 years   ", "1 year   ", "Time lived in place",
                                   "15   ", "9   ", "5   ", "1   ", "None   ", "Experience in local politics", 
                                   "Mostly federal   ", "Majority federal   ", "Majority local   ", "Mostly local   ", "Position on federalism"))

# get values Germany
out[2:5,5:7] <- regDE[26:29,c(6,10:11)]
out[7:11,5:7] <- regDE[21:25,c(6,10:11)]
out[13:18,5:7] <- regDE[15:20,c(6,10:11)]
out[20:22,5:7] <- regDE[12:14,c(6,10:11)]
out[24:28,5:7] <- regDE[7:11,c(6,10:11)]
out[30:31,5:7] <- regDE[5:6,c(6,10:11)]
out[33:36,5:7] <- regDE[1:4,c(6,10:11)]

out[38:41,5:7] <- regDE[55:58,c(6,10:11)]
out[43:47,5:7] <- regDE[50:54,c(6,10:11)]
out[49:54,5:7] <- regDE[44:49,c(6,10:11)]
out[56:58,5:7] <- regDE[41:43,c(6,10:11)]
out[60:64,5:7] <- regDE[36:40,c(6,10:11)]
out[66:67,5:7] <- regDE[34:35,c(6,10:11)]
out[69:72,5:7] <- regDE[30:33,c(6,10:11)]

# get values for the United States
out[74:77,5:7] <- regUS[23:26,c(6,10:11)]
out[79:83,5:7] <- regUS[18:22,c(6,10:11)]
out[85:90,5:7] <- regUS[12:17,c(6,10:11)]
out[92:96,5:7] <- regUS[7:11,c(6,10:11)]
out[98:99,5:7] <- regUS[5:6,c(6,10:11)]
out[101:104,5:7] <- regUS[1:4,c(6,10:11)]

out[106:109,5:7] <- regUS[49:52,c(6,10:11)]
out[111:115,5:7] <- regUS[44:48,c(6,10:11)]
out[117:122,5:7] <- regUS[38:43,c(6,10:11)]
out[124:128,5:7] <- regUS[33:37,c(6,10:11)]
out[130:131,5:7] <- regUS[31:32,c(6,10:11)]
out[133:136,5:7] <- regUS[27:30,c(6,10:11)]

# plot
p <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out[out$country %in% c("Germany", "United States"),], aes(x=attribute, y=b, shape=sample, col=sample), 
             size=2.5, position=position_dodge(width=0.8)) +
  geom_errorbar(data=out[out$country %in% c("Germany", "United States"),], aes(x=attribute, ymin=lb, ymax=ub, col=sample), width=0,
                position=position_dodge(width=0.8), linewidth=1) + 
  scale_color_manual(values = c("burlywood3","burlywood4"), name="") +
  scale_shape_manual(values = c(15,17), name="") +
  scale_y_continuous(breaks = seq(0.3,0.7,0.05)) +
  facet_wrap(~country, scales="free_x") +
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=13,
                                   face = c(rep("plain",4),"bold",
                                            rep("plain",2), "bold",
                                            rep("plain",5), "bold",
                                            rep("plain",3), "bold",
                                            rep("plain",6), "bold",
                                            rep("plain",5), "bold",
                                            rep("plain",4),"bold")),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))
p

ggsave("FigureA4_Attachments_Conjoint_full.png",p,width=10,height=9)





# Figure A5: Conjoint AMCE ------------------------------------------------

### prepare output container
out <- data.frame(matrix(nrow=66, ncol=7))
colnames(out) <- c("country", "sample", "dimension", "attribute", "b", "lb", "ub")
out$country <- factor(rep(1:3,each=22), levels = 1:3,
                      labels = c("Germany", "United States", "Difference in AMCE"))
out$sample <- factor(c(rep(rep(1:2,each=11),2), rep(3:4, each=11)), levels = 1:4,
                     labels = c("Strong", "Weak", "Germany", "United States"))
out$dimension <- factor(rep(c(1,1,1,2,2,2,2,3,3,3,3),6), levels = 1:3,
                        labels = c("Position on federalism", "Experience in local politics", "Time lived in place"))
out$attribute <- factor(rep(11:1,6), levels = 1:11,
                        labels = c("20+ years   ", "10-15 years   ", "1-5 years   ", "Time lived in place",
                                   "9-15 years   ","1-5  years   ", "None   ", "Experience in local politics", 
                                   "Rather local   ","Rather federal   ", "Position on federalism"))

### AMCEs
# Germany
reg <- cj(data = de.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + C_party + 
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "amce",
          by=~loc.attachDICH)

out[2,5:7] <- 0
out[3,5:7] <- reg[44,c(6,10,11)]
out[5,5:7] <- 0
out[6:7,5:7] <- reg[41:42,c(6,10,11)]
out[9,5:7] <- 0
out[10:11,5:7] <- reg[38:39,c(6,10,11)]

out[13,5:7] <- 0
out[14,5:7] <- reg[22,c(6,10,11)]
out[16,5:7] <- 0
out[17:18,5:7] <- reg[19:20,c(6,10,11)]
out[20,5:7] <- 0
out[21:22,5:7] <- reg[16:17,c(6,10,11)]


# United States
reg <- cj(data = us.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + 
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "amce",
          by=~loc.attachDICH)

out[24,5:7] <- 0
out[25,5:7] <- reg[38,c(6,10,11)]
out[27,5:7] <- 0
out[28:29,5:7] <- reg[35:36,c(6,10,11)]
out[31,5:7] <- 0
out[32:33,5:7] <- reg[32:33,c(6,10,11)]

out[35,5:7] <- 0
out[36,5:7] <- reg[19,c(6,10,11)]
out[38,5:7] <- 0
out[39:40,5:7] <- reg[16:17,c(6,10,11)]
out[42,5:7] <- 0
out[43:44,5:7] <- reg[13:14,c(6,10,11)]



### Difference in AMCEs
# Germany
reg <- cj(data = de.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + C_party +
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "amce_difference",
          by=~loc.attachDICH)

out[46,5:7] <- 0
out[47,5:7] <- reg[4,c(7,11,12)]
out[49,5:7] <- 0
out[50:51,5:7] <- reg[6:7,c(7,11,12)]
out[53,5:7] <- 0
out[54:55,5:7] <- reg[14:15,c(7,11,12)]

# United States
reg <- cj(data = us.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + 
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "amce_difference",
          by=~loc.attachDICH)

out[57,5:7] <- 0
out[58,5:7] <- reg[4,c(7,11,12)]
out[60,5:7] <- 0
out[61:62,5:7] <- reg[6:7,c(7,11,12)]
out[64,5:7] <- 0
out[65:66,5:7] <- reg[12:13,c(7,11,12)]


### create new variables for faceting
out$facet_country <- NA
out$facet_country[out$country == "Germany" | out$sample == "Germany"] <- 1
out$facet_country[out$country == "United States" | out$sample == "United States"] <- 2
out$facet_country <- factor(out$facet_country, levels = 1:2,
                            labels = c("Germany", "United States"))

out$facet_analysis <- NA
out$facet_analysis[out$country %in% c("Germany", "United States")] <- 1
out$facet_analysis[out$country == "Difference in AMCE"] <- 2
out$facet_analysis <- factor(out$facet_analysis, levels = 1:2,
                             labels = c("AMCE", "Difference in\nAMCE"))

out$hline <- 0

out$sample2 <- as.numeric(out$sample)
out$sample2[out$sample2 == 4] <- 3
out$sample2 <- factor(out$sample2, levels = 1:3,
                      labels = c("Strong", "Weak", "Strong - Weak"))


### Plot
p <- ggplot() + 
  geom_hline(data=out, aes(yintercept = 0), lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, shape=sample2, col=sample2), 
             size=3, position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, col=sample2), width=0,
                position=position_dodge(width=0.8), linewidth=1) + 
  scale_color_manual(values = c("burlywood4","burlywood3", "gray40"), name="") +
  scale_shape_manual(values = c(15,17,19), name="") +
  facet_grid(facet_country ~ facet_analysis, scales="free_x", space="free_x") +
  coord_flip() + 
  ylab("Average Marginal Component Effect") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=13,
                                   face = c(rep("plain",3),"bold",
                                            rep("plain",3), "bold",
                                            rep("plain",2), "bold")),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))
p
ggsave("FigureA5_Attachments_Conjoint_AMCE.png",p,width=10,height=6)






# Figure A6: correlations with control for place identity -----------------

### Write a function
regfun <- function(dv){
  # define outcome
  de$dv <- as.numeric(de[,dv])
  us$dv <- as.numeric(us[,dv])
  
  # regressions
  outUS1 <- summary(lm_robust(dv ~ loc.attach, 
                              data=us, weights=us$weight))$coef[2,1:2]
  outUS2 <- summary(lm_robust(dv ~ loc.attach + identity_rururb, 
                              data=us, weights=us$weight))$coef[2,1:2]
  
  # convert to SD changes
  outUS1.SD <- sqrt(weighted.var(us$loc.attach, na.rm=T, w=us$weight))*outUS1 / sqrt(weighted.var(us$dv, na.rm=T, w=us$weight))
  outUS2.SD <- sqrt(weighted.var(us$loc.attach, na.rm=T, w=us$weight))*outUS2 / sqrt(weighted.var(us$dv, na.rm=T, w=us$weight))
  
  # return
  return(rbind(outUS1.SD, outUS2.SD))
}


### output container
out <- data.frame(matrix(nrow=42,ncol=6))
colnames(out) <- c("analysis","outcome", "level", "sample", "b", "se")
out$analysis <- factor(rep(c(1,1,1,1,2,2), 7), levels = 1:2,
                       labels = c("level", "Difference federal vs. local"))
out$outcome <- factor(rep(7:1, each=6), levels = 1:7,
                      labels = c("Identification",
                                 "Importance of turnout",
                                 "News consumption",
                                 "Discuss with family",
                                 "Gov't has impact on respondent",
                                 "Respondent can influence gov't",
                                 "Trust in gov't"))
out$level <- factor(rep(rep(1:3,each=2), 7), levels = 1:3,
                    labels = c("Federal", "Local", "Difference"))
out$sample <- factor(rep(1:2), levels = 1:2,
                     labels = c("excl. identity", "incl. identity"))


### Get values
out[1:2,5:6] <- regfun(dv = "govlev_trust_federal")
out[3:4,5:6] <- regfun(dv = "govlev_trust_local")
out[5:6,5:6] <- regfun(dv = "govlev_trustDiff_local")

out[7:8,5:6] <- regfun(dv = "efficacy_federal")
out[9:10,5:6] <- regfun(dv = "efficacy_local")
out[11:12,5:6] <- regfun(dv = "efficacyDiff_local")

out[13:14,5:6] <- regfun(dv = "govlev_effect_federal")
out[15:16,5:6] <- regfun(dv = "govlev_effect_local")
out[17:18,5:6] <- regfun(dv = "govlev_effectDiff_local")

out[19:20,5:6] <- regfun(dv = "discuss_federal")
out[21:22,5:6] <- regfun(dv = "discuss_local")
out[23:24,5:6] <- regfun(dv = "discuss_difference")

out[25:26,5:6] <- regfun(dv = "news_federal")
out[27:28,5:6] <- regfun(dv = "news_local")
out[29:30,5:6] <- regfun(dv = "news_difference")

out[31:32,5:6] <- NA
out[33:34,5:6] <- NA
out[35:36,5:6] <- regfun(dv = "turnout_level")

out[37:38,5:6] <- NA
out[39:40,5:6] <- NA
out[41:42,5:6] <- regfun(dv = "identity_country_local")



### Plot
p1 <- ggplot() +
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out[out$analysis=="level",], aes(x=outcome, y=b, col=level, shape=level),
             position = position_dodge(width=0.6), size=2.5) +
  geom_errorbar(data=out[out$analysis=="level",], aes(x=outcome, ymin=b-1.96*se, ymax=b+1.96*se, col=level),
                position = position_dodge(width=0.6), size=1, width=0.25) +
  scale_y_continuous(limits = c(-0.15,0.4), breaks = seq(-0.1,0.4,0.1)) +  
  facet_wrap(~sample) +
  scale_color_manual(values=c("burlywood3", "burlywood4"), name="") +
  scale_shape_manual(values=c(19,17), name="") +
  ylab("Estimated relationship with local attachments\n(standard deviations)") +
  coord_flip() +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        #axis.title.x = element_text(size=14, face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size=12), 
        axis.text = element_text(size=12))  

p2 <- ggplot() +
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out[out$analysis=="Difference federal vs. local",], aes(x=outcome, y=b, col=sample, shape=sample),
             position = position_dodge(width=0.6), size=2.5) +
  geom_errorbar(data=out[out$analysis=="Difference federal vs. local",], aes(x=outcome, ymin=b-1.96*se, ymax=b+1.96*se, col=sample),
                position = position_dodge(width=0.6), size=1, width=0.25) +
  scale_y_continuous(limits = c(-0.25,0.1)) +
  facet_wrap(~analysis) +
  scale_color_manual(values=c("chocolate2", "chocolate4"), name="") +
  scale_shape_manual(values=c(19,17), name="") +
  ylab("Estimated relationship with local attachments\n(standard deviations)") +
  coord_flip() +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        #axis.title.x = element_text(size=14, face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size=12), 
        axis.text.x = element_text(size=12),
        axis.text.y = element_blank())  

p <- ggarrange(p1,p2,widths=c(0.67,0.33))
p

ggsave("FigureA6_Attachments_Differences_URidentity.png",p,width=10,height=6)





# Figure A7: conjoint results by local identity ---------------------------

### combine local attachments and local identity scale (US only)
us.conjoint$identity_combined <- NA
us.conjoint$identity_combined[us.conjoint$loc.attachDICH == "weak" & us.conjoint$identity_rururb_partDICH == 0] <- 1
us.conjoint$identity_combined[us.conjoint$loc.attachDICH == "weak" & us.conjoint$identity_rururb_partDICH == 1] <- 2
us.conjoint$identity_combined[us.conjoint$loc.attachDICH == "strong" & us.conjoint$identity_rururb_partDICH == 0] <- 3
us.conjoint$identity_combined[us.conjoint$loc.attachDICH == "strong" & us.conjoint$identity_rururb_partDICH == 1] <- 4
us.conjoint$identity_combined <- factor(us.conjoint$identity_combined, levels = 1:4,
                                        labels = c("weak attach, weak ident",
                                                   "weak attach, strong ident",
                                                   "strong attach, weak ident",
                                                   "strong attach, strong ident"))


### prepare output container
out <- data.frame(matrix(nrow=88+22, ncol=8))
colnames(out) <- c("statistic","attachments", "identity", "dimension", "attribute", "b", "lb", "ub")
out$statistic <- factor(c(rep(1:2,each=44), rep(3,22)), levels = 1:3,
                        labels = c("Marginal Means", "AMCE", "Difference in\nAMCEs"))
out$attachments <- factor(c(rep(rep(1:2,each=22),2),rep(3,22)), levels = 1:3,
                          labels = c("Weak attachments", "Strong attachments",
                                     "Strong - Weak"))
out$identity <- factor(c(rep(rep(1:2,each=11),4), rep(1:2,each=11)), levels = 1:2,
                       labels = c("Weak\nrural/urban identity", "Strong\nrural/urban identity"))
out$dimension <- factor(rep(c(1,1,1,2,2,2,2,3,3,3,3),10), levels = 1:3,
                        labels = c("Position on federalism", "Experience in local politics", "Time lived in place"))
out$attribute <- factor(rep(11:1,10), levels = 1:11,
                        labels = c("20+ years   ", "10-15 years   ", "1-5 years   ", "Time lived in place",
                                   "9-15 years   ","1-5  years   ", "None   ", "Experience in local politics", 
                                   "Rather local   ","Rather federal   ", "Position on federalism"))


### Marginal means
# compute
reg <- cj(data = us.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + 
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "mm",
          by=~identity_combined)

# weak attachments, weak identity
out[2:3,6:8] <- reg[18:19, c(6,10,11)]
out[5:7,6:8] <- reg[15:17, c(6,10,11)]
out[9:11,6:8] <- reg[12:14, c(6,10,11)]

# weak attachments, strong identity
out[13:14,6:8] <- reg[37:38, c(6,10,11)]
out[16:18,6:8] <- reg[34:36, c(6,10,11)]
out[20:22,6:8] <- reg[31:33, c(6,10,11)]

# strong attachments, weak identity
out[24:25,6:8] <- reg[56:57, c(6,10,11)]
out[27:29,6:8] <- reg[53:55, c(6,10,11)]
out[31:33,6:8] <- reg[50:52, c(6,10,11)]

# strong attachments, strong identity
out[35:36,6:8] <- reg[75:76, c(6,10,11)]
out[38:40,6:8] <- reg[72:74, c(6,10,11)]
out[42:44,6:8] <- reg[69:71, c(6,10,11)]


### AMCEs
# compute
reg <- cj(data = us.conjoint, formula = C_choice ~ C_age + C_gender + C_occup + 
            C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
          estimate = "amce",
          by=~identity_combined)

# weak attachments, weak identity
out[46:47,6:8] <- reg[18:19, c(6,10,11)]
out[49:51,6:8] <- reg[15:17, c(6,10,11)]
out[53:55,6:8] <- reg[12:14, c(6,10,11)]

# weak attachments, strong identity
out[57:58,6:8] <- reg[37:38, c(6,10,11)]
out[60:62,6:8] <- reg[34:36, c(6,10,11)]
out[64:66,6:8] <- reg[31:33, c(6,10,11)]

# strong attachments, weak identity
out[68:69,6:8] <- reg[56:57, c(6,10,11)]
out[71:73,6:8] <- reg[53:55, c(6,10,11)]
out[75:77,6:8] <- reg[50:52, c(6,10,11)]

# strong attachments, strong identity
out[79:80,6:8] <- reg[75:76, c(6,10,11)]
out[82:84,6:8] <- reg[72:74, c(6,10,11)]
out[86:88,6:8] <- reg[69:71, c(6,10,11)]


### Diff. in AMCE
# compute
reg1 <- cj(data = us.conjoint[us.conjoint$identity_combined %in% c("weak attach, weak ident", "strong attach, weak ident"),], 
           formula = C_choice ~ C_age + C_gender + C_occup + 
             C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
           estimate = "amce_difference",
           by=~identity_combined)
reg2 <- cj(data = us.conjoint[us.conjoint$identity_combined %in% c("weak attach, strong ident", "strong attach, strong ident"),], 
           formula = C_choice ~ C_age + C_gender + C_occup + 
             C_timelived_combined + C_localexp_combined + C_federalism_combined, id = ~ResponseId,
           estimate = "amce_difference",
           by=~identity_combined)

# weak identity
out[90,6] <- 0
out[91,6:8] <- reg1[4,c(7,11,12)]
out[93,6] <- 0
out[94:95,6:8] <- reg1[6:7,c(7,11,12)]
out[97,6] <- 0
out[98:99,6:8] <- reg1[12:13,c(7,11,12)]

# strong identity
out[101,6] <- 0
out[102,6:8] <- reg2[4,c(7,11,12)]
out[104,6] <- 0
out[105:106,6:8] <- reg2[6:7,c(7,11,12)]
out[108,6] <- 0
out[109:110,6:8] <- reg2[12:13,c(7,11,12)]



### define hline
out$hline <- c(rep(c(0.5,0), each=44), rep(0,22))


### Plot
p <- ggplot() + 
  geom_hline(data=out, aes(yintercept = hline), lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, shape=attachments, col=attachments), 
             size=3, position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, col=attachments), width=0,
                position=position_dodge(width=0.8), linewidth=1) + 
  scale_color_manual(values = c("burlywood3","burlywood4", "gray40"), name="") +
  scale_shape_manual(values = c(15,17,19), name="") +
  facet_grid(identity~statistic, scales="free_x") +
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=13,
                                   face = c(rep("plain",3),"bold",
                                            rep("plain",3), "bold",
                                            rep("plain",2), "bold")),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))
p

ggsave("FigureA7_Attachments_Conjoint_identityUS.png",p,width=11,height=6)






# Figure A8: open-ended place attachment question -------------------------

### Output container
out <- data.frame(matrix(nrow=9, ncol=4))
colnames(out) <- c("class", "variable", "n", "share")
out$class <- factor(c(0, rep(1,3), 2, rep(3,4)), levels = 0:3,
                    labels = c("Any positive answer", "Social ties", "Biographical ties", "Physical-cultural ties"))
out$variable <- factor(1:9, levels = 1:9,
                       labels = c("Any\npositive answer",
                                  "Social ties\nall", "Family", "Friends/\nAcquaintances",
                                  "Biographical\nties all",
                                  "Physical-cultural\nties all", "Physical", "Cultural", "Natural"))

### Get values
# any answer given
out[1,3] <- table(I(de$text_place_answer == "YES"))[2]
out[1,4] <- prop.table(table(I(de$text_place_answer == "YES")))[2]

# social ties combined
de$text_place_socialties <- ifelse(de$text_place_weakties == 1 | de$text_place_strongties == 1, 1, 0)
out[2,3] <- table(de$text_place_socialties)[2]
out[2,4] <- prop.table(table(de$text_place_socialties))[2]

# strong ties
out[3,3] <- table(de$text_place_strongties)[2]
out[3,4] <- prop.table(table(de$text_place_strongties))[2]

# weak ties
out[4,3] <- table(de$text_place_weakties)[2]
out[4,4] <- prop.table(table(de$text_place_weakties))[2]

# biographical ties
out[5,3] <- table(de$text_place_biographical)[2]
out[5,4] <- prop.table(table(de$text_place_biographical))[2]

# physical-cultural ties combined
de$text_place_physculturalties <- ifelse(de$text_place_physical == 1 | de$text_place_cultural == 1 | de$text_place_natural == 1, 1, 0)
out[6,3] <- table(de$text_place_physculturalties)[2]
out[6,4] <- prop.table(table(de$text_place_physculturalties))[2]

# physical
out[7,3] <- table(de$text_place_physical)[2]
out[7,4] <- prop.table(table(de$text_place_physical))[2]

# cultural
out[8,3] <- table(de$text_place_cultural)[2]
out[8,4] <- prop.table(table(de$text_place_cultural))[2]

# natural
out[9,3] <- table(de$text_place_natural)[2]
out[9,4] <- prop.table(table(de$text_place_natural))[2]


### define alpha
out$alpha <- factor(c(1,1,0,0,1,1,0,0,0), levels = c(0,1),
                    labels = c("component", "all"))


### plot
p <- ggplot() +
  geom_bar(data=out[out$class != "Any positive answer",], 
           aes(x=variable, y=share*100, fill=class, alpha=alpha),
           position=position_dodge(), stat="identity") +
  scale_alpha_manual(values=c(0.6,1), name="") +
  scale_fill_manual(values = c("burlywood4", "burlywood3", "burlywood2"), name="") +
  scale_y_continuous(breaks=seq(0,100,10)) +
  ylab("Share of respondents (%)") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=12, face="bold"),
        axis.text.x = element_text(size=12, face=c("bold", "plain", "plain",
                                                   "bold",
                                                   "bold", "plain", "plain", "plain")),
        axis.text.y = element_text(size=12),
        legend.position = "none",
        legend.text = element_text(size=12))  
p
ggsave("FigureA8_OpenEnds_Q1.png",p,width=10.5,height=5)





# Figure A9: open-ended place identity question ---------------------------

### output container
out <- data.frame(matrix(nrow=10, ncol=6))
colnames(out) <- c("variable", "label", "N_low", "N_high", "Share_low", "Share_high")
out$variable <- factor(rep(1:5,each=2), levels = 1:5,
                       labels = c("Any commonality", "Identity-narrow", "Identity-broad", "Identity-combined", "Non-identity"))
out$label <- factor(rep(1:2,5), levels=1:2,
                    labels = c("no", "yes"))

### get values
# answer given?
out[1:2,3:4] <- table(de$text_identity_answer, de$loc.attachDICH)
out[1:2,5:6] <- prop.table(table(de$text_identity_answer, de$loc.attachDICH),2)

# identity narrow
out[3:4,3:4] <- table(de$text_identity_narrow, de$loc.attachDICH)
out[3:4,5:6] <- prop.table(table(de$text_identity_narrow, de$loc.attachDICH),2)

# identity broad
out[5:6,3:4] <- table(de$text_identity_broad, de$loc.attachDICH)
out[5:6,5:6] <- prop.table(table(de$text_identity_broad, de$loc.attachDICH),2)

# identity combined
de$text_identity_broad_combined <- ifelse(de$text_identity_narrow == 1 | de$text_identity_broad == 1, 1, 0)
out[7:8,3:4] <- table(de$text_identity_broad_combined, de$loc.attachDICH)
out[7:8,5:6] <- prop.table(table(de$text_identity_broad_combined, de$loc.attachDICH),2)

# non-identity
out[9:10,3:4] <- table(de$text_identity_nonidentity, de$loc.attachDICH)
out[9:10,5:6] <- prop.table(table(de$text_identity_nonidentity, de$loc.attachDICH),2)


### output table
print(xtable(out, caption="Responses to identity question, by strength of place attachment"), 
      include.rownames=FALSE, caption.placement = "top")



### turn the table into a plot
# reshape data
toplot1 <- subset(out, select = c(variable, label, Share_low), label == "yes")
colnames(toplot1)[3] <- "share"
toplot2 <- subset(out, select = c(variable, label, Share_high), label == "yes")
colnames(toplot2)[3] <- "share"
toplot1$loc.attachDICH <- 1
toplot2$loc.attachDICH <- 2
toplot <- rbind(toplot1,toplot2)
toplot$loc.attachDICH <- factor(toplot$loc.attachDICH, levels = 1:2,
                     labels = c("Weak attachments", "Strong attachments"))

# plot
p <- ggplot() + 
  geom_bar(data=toplot, aes(x=variable, y=share*100, fill=loc.attachDICH),
           stat="identity", position=position_dodge()) +
  scale_x_discrete(limits = rev(levels(toplot$variable))) +
  scale_y_continuous(breaks=seq(0,80,10)) +
  scale_fill_manual(values = c("burlywood3","burlywood4"), name="") +
  coord_flip() +
  xlab("") +
  ylab("Share of respondents (%)") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title = element_text(size=12, face="bold"),
        legend.position = "bottom",
        legend.text = element_text(size=12), 
        axis.text = element_text(size=12))
p
ggsave("FigureA9_OpenEnds_Q2.png",p,width=7,height=4)





# Table A4: regional US sample details ---------------------------

### unweighted
app %>%
  mutate(age_30_39 = age_cat_i=="age_30_39") %>%
  mutate(age_70p = age_cat_i=="age_70p") %>%
  mutate(b = race_2=="Black or African American") %>%
  mutate(w = race_1=="White") %>%
  mutate(male = gender=="M") %>%
  mutate(IN = state_name=="Indiana") %>%
  mutate(OH = state_name=="Ohio") %>%
  mutate(PA = state_name=="Pennsylvania") %>%
  mutate(WV = state_name=="West Virginia") %>%
  mutate(KY = state_name=="Kentucky") %>%
  mutate(dem = pid3=="D") %>%
  mutate(rep = pid3=="R") %>%
  mutate(ed = educ<5) %>%
  dplyr::summarize(across(c(counties, age_30_39:ed), ~mean(., na.rm=T)))

### weighted
app %>%
  mutate(age_30_39 = age_cat_i=="age_30_39") %>%
  mutate(age_70p = age_cat_i=="age_70p") %>%
  mutate(b = race_2=="Black or African American") %>%
  mutate(w = race_1=="White") %>%
  mutate(male = gender=="M") %>%
  mutate(IN = state_name=="Indiana") %>%
  mutate(OH = state_name=="Ohio") %>%
  mutate(PA = state_name=="Pennsylvania") %>%
  mutate(WV = state_name=="West Virginia") %>%
  mutate(KY = state_name=="Kentucky") %>%
  mutate(dem = pid3=="D") %>%
  mutate(rep = pid3=="R") %>%
  mutate(ed = educ<5) %>%
  dplyr::summarize(across(c(counties, age_30_39:ed), ~weighted.mean(., w=weight, na.rm=T)))









